﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using Meta.Numerics.Statistics.Distributions;

namespace ExcelAddIn1
{
    /**
     * Implementation of the Kruskal-Wallis test.
     */
    class Test_KruskalWallis : Test
    {
        private class RankItem : IComparable
        {
            public double value;
            public String dataset;
            public double rank;
        
            public RankItem(double value, String dataset)
            {
                this.value = value;
                this.dataset = dataset;
                this.rank = 0;
            }
        
            public int CompareTo(Object obj)
            {
                RankItem ri = (RankItem)obj;
                if (value < ri.value)
                {
                    return -1;
                }
                else if (value > ri.value)
                {
                    return 1;
                }
                else
                {
                    return 0;
                }
            }
        }

        private double[,] R005 = {
            {6.085, 8.331, 9.798, 10.881, 11.734, 12.435, 13.027, 13.538, 13.988},
            {4.501, 5.910, 6.825, 7.502, 8.037, 8.478, 8.852, 9.177, 9.462},
            {3.927, 5.040, 5.757, 6.287, 6.707, 7.053, 7.347, 7.602, 7.826},
            {3.635, 4.602, 5.219, 5.673, 6.033, 6.330, 6.582, 6.801, 6.995},
            {3.461, 4.339, 4.896, 5.305, 5.629, 5.895, 6.122, 6.319, 6.493},
            {3.344, 4.165, 4.681, 5.060, 5.359, 5.606, 5.815, 5.998, 6.158},
            {3.261, 4.041, 4.529, 4.886, 5.167, 5.399, 5.596, 5.767, 5.918},
            {3.199, 3.949, 4.415, 4.755, 5.024, 5.244, 5.432, 5.595, 5.738},
            {3.151, 3.877, 4.327, 4.654, 4.912, 5.124, 5.304, 5.461, 5.598},
            {3.113, 3.820, 4.256, 4.574, 4.823, 5.028, 5.202, 5.353, 5.486},
            {3.081, 3.773, 4.199, 4.508, 4.748, 4.947, 5.116, 5.263, 5.395},
            {3.055, 3.734, 4.151, 4.453, 4.690, 4.884, 5.049, 5.192, 5.318},
            {3.033, 3.701, 4.111, 4.407, 4.639, 4.829, 4.990, 5.130, 5.253},
            {3.014, 3.673, 4.076, 4.367, 4.595, 4.782, 4.940, 5.077, 5.198},
            {2.998, 3.649, 4.046, 4.333, 4.557, 4.741, 4.896, 5.031, 5.150},
            {2.984, 3.628, 4.020, 4.303, 4.524, 4.705, 4.858, 4.991, 5.108},
            {2.971, 3.609, 3.997, 4.276, 4.494, 4.673, 4.824, 4.955, 5.071},
            {2.960, 3.593, 3.977, 4.253, 4.469, 4.645, 4.794, 4.924, 5.038},
            {2.950, 3.578, 3.958, 4.232, 4.445, 4.620, 4.768, 4.895, 5.008},
            {2.943, 3.566, 3.943, 4.214, 4.425, 4.599, 4.745, 4.871, 4.982},
            {2.935, 3.554, 3.928, 4.197, 4.407, 4.578, 4.723, 4.848, 4.958},
            {2.927, 3.543, 3.915, 4.182, 4.389, 4.559, 4.703, 4.827, 4.936},
            {2.920, 3.533, 3.902, 4.167, 4.374, 4.542, 4.685, 4.808, 4.916}
        };

        public override bool EnableAlphaField()
        {
            return false;
        }

        public Test_KruskalWallis(DataContainer data)
        {
            this.data = data;
            res = new List<String>();
        }

        public override string GetInfo()
        {
            String s = "Performs a Kruskal-Wallis test to check for significant differences between three or more datasets. ";
            s += "The test does not require normally distributed data.";
            return s;
        }

        public override void RunTest()
        {
            res.Add("Kruskal-Wallis test");

            //Error check
            if (data.GetNoSets() < 3)
            {
                res.Add("At least three datasets is required!");
                return;
            }

            //Get data sets
            int totalN = 0;
            for (int s = 0; s < data.GetNoSets(); s++)
            {
                totalN += data.GetDataSet(s).GetN();
            }

            //Create the rankings table
            List<RankItem> rankTab = new List<RankItem>(totalN);
            //Fill table with all datasets
            for (int s = 0; s < data.GetNoSets(); s++)
            {
                DataSet d = data.GetDataSet(s);
                for (int i = 0; i < d.GetN(); i++)
                {
                    rankTab.Add(new RankItem(d.GetValue(i), d.GetName()));
                }
            }

            //Sort the table
            rankTab.Sort();
            //Fill in ranks
            for (int i = 0; i < totalN; i++)
            {
                rankTab[i].rank = i + 1;
            }
            //Fix rank duplicates
            double prev = rankTab[0].value;
            double sumRank = rankTab[0].rank;
            int cnt = 1;
            for (int i = 1; i < totalN; i++)
            {
                double iVal = rankTab[i].value;
                if (iVal == prev)
                {
                    //Equal values. Find a mean rank.
                    sumRank += rankTab[i].rank;
                    cnt++;
                }
                else
                {
                    if (cnt > 1)
                    {
                        //Duplicates found. Update rank to a mean rank.
                        double meanR = sumRank / (double)(cnt);
                        for (int d = i - 1; d >= i - cnt; d--)
                        {
                            rankTab[d].rank = meanR;
                        }
                    }
                    cnt = 1;
                    sumRank = rankTab[i].rank;
                }
                prev = iVal;
            }
            //Rankings table done!

            //Calculate R-scores
            double[] R = new double[data.GetNoSets()];
            for (int s = 0; s < data.GetNoSets(); s++)
            {
                DataSet d = data.GetDataSet(s);
                for (int n = 0; n < totalN; n++)
                {
                    RankItem ri = rankTab[n];
                    if (ri.dataset.Equals(d.GetName()))
                    {
                        R[s] += ri.rank;
                    }
                }
            }

            //Calculate H-score
            double Rsum = 0;
            for (int s = 0; s < data.GetNoSets(); s++)
            {
                DataSet d = data.GetDataSet(s);
                double ns = (double)d.GetN();
                Rsum += Math.Pow(R[s], 2) / ns;
            }
            double H = 12.0 / ((double)totalN * (double)(totalN + 1)) * Rsum - 3.0 * (double)(totalN + 1);

            //Find Critical H
            int DF = data.GetNoSets() - 1;
            ChiSquaredDistribution csd = new ChiSquaredDistribution(DF);
            double Hc = csd.InverseRightProbability(alpha);

            //Results
            res.Add("DoF;" + DF + ";;-");
            res.Add("α;" + alpha);
            res.Add("H-score;" + H.ToString("F2"));
            res.Add("H-crit;" + Hc.ToString("F2")); 
            if (H > Hc)
            {
                res.Add("Result;Significant difference");
            }
            else
            {
                res.Add("Result;No difference");
            }
            res.Add(";;;-");

            //Perform Tukey's post-test
            res.Add("Tukey's post-test");
            int k = data.GetNoSets();

            //Calculate total sum of squares and total n.
            double ntot = 0;
            double totSum = 0.0;
            double totSqSum = 0.0;
            for (int s = 0; s < k; s++)
            {
                DataSet d = data.GetDataSet(s);
                ntot += (double)d.GetN();
                totSum += d.GetSum();
                totSqSum += d.GetSquaredSum();
            }
            double SST = totSqSum - Math.Pow(totSum, 2) / ntot;

            //Calculate grand mean
            double grandMean = 0.0;
            for (int s = 0; s < k; s++)
            {
                DataSet d = data.GetDataSet(s);
                grandMean += d.GetMean() * (double)d.GetN();
            }
            grandMean = grandMean / ntot;

            //Calculate SSB and SSW
            double SSB = 0.0;
            for (int s = 0; s < k; s++)
            {
                DataSet d = data.GetDataSet(s);
                SSB += (double)d.GetN() * Math.Pow(d.GetMean() - grandMean, 2);
            }
            double SSW = SST - SSB;

            //Calculate F-score
            double MSB = SSB / (double)(k - 1);
            double MSW = SSW / (ntot - (double)k);
            double F = MSB / MSW;

            //Get Critical R
            int DF1 = k;
            int DF2 = (k - 1) * (data.GetDataSet(0).GetN() - 1);
            if (DF1 > 24 || DF2 > 10)
            {
                res.Add("Tukey's post-test does not support DoF over 24:10!");
                return;
            }
            double q = R005[DF2 - 2, DF1 - 2];
            double Rc = q * Math.Sqrt(MSW / (double)data.GetDataSet(0).GetN());

            res.Add("R-crit;" + Rc.ToString("F2") + ";;-");
            for (int i = 0; i < data.GetNoSets(); i++)
            {
                for (int j = 0; j < data.GetNoSets(); j++)
                {
                    if (i != j && i < j)
                    {
                        DataSet d1 = data.GetDataSet(i);
                        DataSet d2 = data.GetDataSet(j);
                        double absD = Math.Abs(d1.GetMean() - d2.GetMean());

                        if (absD > Rc)
                        {
                            res.Add("Rs(" + (i + 1) + "," + (j + 1) + ");" + absD.ToString("F2") + ";Difference");
                        }
                        else
                        {
                            res.Add("Rs(" + (i + 1) + "," + (j + 1) + ");" + absD.ToString("F2") + ";No difference");
                        }
                    }
                }
            }
            res.Add(";;;-");
        }
    }
}
